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' Abstract. We initiate numerical studies of differentially rotating magnetised (proto) neutron stars by studying - 

through construction from first principles - the coupling between an assumed differential rotation and an impressed 

, magnetic field. For a perfect incompressible, homogeneous, non-dissipative fiuid sphere immersed in an ambient 

plasma, we solve the (coupled) azimuthal components of the Navier-Stokes equation and the Maxwell induction 
equation. The assumed time-independent poloidal field lines get dragged by the rotating fluid and produce toroidal 
magnetic fields. Surface magnetic fields take away energy redistributing the angular momentum to produce rigid 
rotation along poloidal field lines. Due to absence of viscous dissipation, sustained torsional oscillations are set 
up within the star. However, the perpetual oscillations of neighbouring 'closed' field lines get increasingly out of 
phase with time, leading to structure build up as in Liu & Shapiro (2004) implying the importance of taking into 
account diffusion (Spruit 1999) for realistic modeling. 

' Key words, gravitational waves - MHD - starsmeutron - stars:rotation - stars: magnetic fields 

' 1. Introduction 

■ Absence of short period (< 1 ms) newborn radio pulsars (see e.g. Bhattacharya & Van den Heuvel 1991) is a long 
' standing conundrum of pulsar physics. Mechanisms for spinning down the protoneutron star - the object left behind 
I immediately (3-4 s) after a supernova explosion or a binary neutron star merger event - over dynamical timescales, 
' therefore, assume significance in modeling these objects. The proposed mechanisms of spin-down are through the 
. emission of gravitational waves (e.g. Lindblom, Owen & Morsink 1998), or through magnetic braking or viscous 

damping (Heger et al 2003). Other processes that rely on instabilities (e.g. Andersson 2003; Spruit & Phinney 1998) 
. for spinning down the star have also received attention. An important ingredient in all these mechanisms is the 

■ differential rotation of the object. Differential rotation is expected to be strong in the first few seconds after the 
supernova collapse (Akiyama et al. 2003) - this can be the case even following a neutron star-neutron star binary 
merger (Zwerger & MuUer 1997; Rampp, MuUer & Ruffert 1998), although such a binary merger results in a black- 
hole formation, following a short lived "hypermassive" neutron star stage (e.g. Shibata & Uryu 2000, Rasio & Shapiro 
1999). Amongst the above mechanisms, the spin down of neutron stars through gravitational wave emission as a 
consequence of r-mode oscillations have aroused particular interest in recent years (Lindblom, Owen & Morsink 1998). 
Particularly, there have been suggestions that neutron stars with superfiuid cores and/or solid crusts could be sources 
for gravitational waves, spinning down young neutron stars to the observed pulsar periods (Andersson & Kokkotas 
2001). The interest in r-mode instability is due to the possibility of the frequency and amplitude of the gravitational 
waves being in the observable range of upcoming and operational gravitational wave detectors like LIGO, VIRGO, 
TAMA, GEO (Heger et al 1998). However, r-modes can suffer suppression or damping in the presence of hyperon bulk 
viscosity (Jones 2001; there nonetheless exist "windows of opportunity" - Andersson & Kokkotas 2001) or magnetic 
fields (RezzoUa, Lamb & Shapiro 2000; RezzoUa et al 2001a,b), in addition, it was recently shown (Arras et al. 2003) 
that the saturation energies of r-modes are extremely small to account for the spin-down to periods of the order of 
10-100 ms. The damping timescales of differential rotation, varies from 10"'^ s for r-modes, to 10^ s for Alfven effects 
or 10^ s for molecular viscosity. Thus, it is of pertinence for pulsar physics, as well as gravitational wave physics, to 
establish which of the above effects dominate and at what timescales they take effect. 



Send offprint requests to: A.V. Thampan 



2 



A.V. Thampan: Differentially rotating magnetised neutron stars: production of toroidal magnetic fields 



Although a magnetic field of 10^^ G is expected by flux conservation of the magnetic field of the progenitor star, the 
exact evolution to a cold rigidly rotating object with a predominantly dipole magnetic field of the above magnitude is 
still a controversial issue (e.g. Chanmugam, Rajasekhar & Young 1995 and references therein). The issue is complicated 
due to observations (for at least one pulsar) of higher mulipole components (Gil & Mitra 2001). The subsequent 
evolution of the magnetic field of the radio pulsar (Srinivasan et al 1990, Ruderman, M. 1991 a,b,c, Blandford, 
Applegate & Hernquist, L. 1983, Sang & Chanmugam 1987, Geppert & Urpin 1994, Urpin & Geppert 1995,1996, 
Konar & Bhattacharya 1997, Konar & Bhattacharya 1999, Mitra, Konar & Bhattacharya 1999, see Bhattacharya 
2002 for a recent review), thus, crucially depends on the magnetic field (and the rotation rate) that the pulsar 
possesses as soon as it is 'switched on' (Vivekanand & Narayan 1981). These effects acquire additional significance 
when gravitational effects are taken into account (RezzoUa, Ahmedov & Miller 2001; Zanotti & RezzoUa 2002). 

A possible mechanism to slow down the pulsar is through the production of interior toroidal magnetic field com- 
ponents (from poloidal components through differential rotation) and back reaction of this on the rotational flow. 
Moffatt (1978) discussed the distortion of an initially uniform magnetic field by differential rotation and showed that 
a net toroidal magnetic field flux over the entire space covering the interior of the star to inflnity, can develop only 
as a result of diffusion and only if the term p(Bp.V)n (where p is the cylindrical radial coordinate. Bp the poloidal 
magnetic field, V the differential operator and the rotation rate) is symmetric about the mid plane. Following such 
an argument, an analytical expression was obtained for the final toroidal field at the asymptotic radial limit for a 
given constant of diffusion. Mestel (1999) explained the need for toroidal currents to maintain the poloidal fields and 
vice versa in the interiors of stars and presented a wave equation for the rotation rate. Spruit (1999) discussed the 
interplay between differential rotation and magnetic fields in stellar interiors, focusing on the inherent instabilities 
that may occur in such systems. 

These works have largely been based on analytical results or analytical extensions of the kinematic dynamo equa- 
tions. There is a need to confirm and extend these results using realistic magneto-hydrodynamic simulations. However, 
such numerical computations have been difficult due to the inherent non-deterministic nature of the problem and due 
to the several numerical instabilities that can occur during the simulation. Only recently have reliable and stable 
numerical schemes been formulated and initial steps towards realistic simulations of these systems been undertaken. 
Shapiro (2000) solved the incompressible, non-turbulent, non-diffusive hydromagnetic dynamo equations both ana- 
lytically and numerically and showed that for assumed seed magnetic fields, viscosity effectively brakes differential 
rotation. It was also shown that an infinite cylindrical star, with an impressed internal monopolar field, can lose 
angular momentum to an external plasma medium, even in the absence of viscosity. The effects of compressibility 
was discussed by Cook, Shapiro & Stephens (2003) who showed that radial flow and shocks become important in the 
evolution. Liu & Shapiro (2004) worked out these effects for relativistic spherical stars with quadrupolar and dipolar 
internal fields. 

Considering, the complexities invloved in these simulations and the high potential for these techniques to be used 
to understand the dynamic behavior of neutron stars and other astrophysical systems, it is prudent to develop different 
numerical schemes to cross verify the results obtained. Moreover, having different reliable schemes gives the flexibility 
to choose the one which can be most easily adapted for a particular application. By far, the best advantage is of 
obtaining different physical insight when the problem is solved for different settings. 

In this work, we present and implement a numerical magneto-hydrodynamic scheme to study the dynamics of 
differentially rotating proto-neutron stars. Assuming the proto-neutron star to be composed of an incompressible 
homogeneous fluid, we evolve the toroidal components of the magnetic field and flow for fixed poloidal components. The 
star is assumed to be spherical in shape and the kinetic energy of rotation is assumed to be far lesser than the magnetic 
energy which itself is far lesser than the gravitational energy. The scheme is tested against numerical dissipation of 
physically conserved quantities and verified by comparing the results obtained with those from similar numerical 
computations (Liu & Shapiro 2004). Although the model used to represent the proto-neutron star is simplistic, we 
show that qualitative insights on the dynamic behavior of such systems can be obtained. The motivation here is to 
present a reliable numerical scheme which may be further developed to simulate realistic neutron stars and other 
astrophysical systems. 

In the next section |j2Jl the details of the technique are presented. In section© the results of the simulations are 
given which are discussed in section Q. Appendix A provides the details of numerical implementation of the equations. 

2. Methodology 

We assume the differentially rotating star to be made up of a homogeneous, incompressible, perfectly conducting fluid. 
The relevant Magnetohydrodynamic (MHD) equations are: 



V • V = 



(1) 
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— + (v • V)v = — VP - V$ + i^V^v + -^^ 

at e 4iT€ 

(2) 

V^* = AwGe (3) 
V • B = (4) 

— = V X (V X B) (5) 

where v is the fluid velocity, P the material pressure, e the matter density, $ the gravitational potential, u the 
kinematic viscosity and B the magnetic field. 

For our purpose here, we assume the star to be spherically symmetric (with a radius R) and immersed in an 
ambient plasma. The initial differential rotation is taken to be constant along cylinders, and the poloidal component 
of the internal magnetic field to be time-independent and given to be a pure dipole: 



1^ {2z'-p') 

i?3 (p2 + ^2)5/2 



where ^ is the magnetic moment. We make use of cylindrical coordinates: p, z (scaled with respect to R) and tp here 
(the spherical symmetry of the star not withstanding, the quantities that are being considered in this problem possess 
axial symmetry in general). Such a magnetic field being divergence and curl free, nonetheless, represents a field in 
vacuum (rather than one that resides inside a medium as is considered here) and contains a singularity at the centre. 
Since this field is time independent, and wc assume a purely axial flow, the central singularity in the magnetic fleld 
does not induce irregular behaviour in the velocity field (even though in the final state, the rotation rate will have a 
singularity at the centre) or the azimutlial magnetic field. In the event of performing more realistic calculations, with 
more realistic magnetic field structures (such as the output of gravitational collapse codes), the numerical techniques 
developed here for treating such a rc;strictive magnetic field strucutre (the central singularity turns up as a stringent 
Courant condition) is expected to prove handy. 

In order to compare with previous calculations in the literature, we use the same initial rotational and toroidal 
magnetic field profiles as that of Shapiro (2000), viz.: 

f2(0,p,z) = lr!o[l + cos(V)] (7) 
5^(0, p, 2) = (8) 

(9) 

where CIq represents a characteristic rotation rate at the centre. 

For our assumptions to be self consistent, we require T A4 <^ W, where T is the rotational kinetic energy, A4, 
the total magnetic energy and W the gravitational potential energy. In other words, for the magnetic field structure 

that we assume here, at time t > 

{BPf + {B'f » {B^f ~ 47rep2jl2 

Since there is no mass in flow to the system and since the system is homogeneous the poloidal equations decouple 
from the toroidal component of the Navier-Stokes equation (Shapiro & Teukolsky 1983). We also assume that viscous 
dissipation, turbulence, convection and diffusion are absent in the system. Another tacit assumption is that the system 
is stable against the innumerable plasma/MHD instabilities (Spruit 1999). 

2.1. Flow and Magnetic Field Parameters 
2.1.1. Formalism 

For our treatment here, we make use of the natural coordinates available to us, viz., one that follows the magnetic 
field lines. For the sake of comparison with the equivalent 1-D problem, we use a modified version of magnetic flux 
coordinates by defining the magnetic field as 



(10) 
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To be consistent with our assumptions of the time independent poloidal field dominating over the time dependent 
toroidal field, the equality above is only approximate. With such a definition, we further choose 

^ ^ ^ (12) 
p ^ p (13) 

to be the first, second and third coordinate respectively. These choices are such that -0 is constant along a magnetic 
field line. The inverse transformations are: 

(14) 

= (15) 



P 




(16) 



In the following paragraphs, we drop the 'carets' from the variables and tacitly understand the variables to be 
dimensionless. With this, treatment, we have for the new coordinate system 



B 


= (^'^,0,5") 




= 




- ^^l/ipB") 


hip 


= = p 


hp 


^ V97p = B/BP 




= ByipiBPf) 




= = 9pv 
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~ 9'P'p{9ipip9pp ~ 9ipp 



\/{BPf 



(17) 
(18) 
(19) 
(20) 
(21) 
(22) 
(23) 
(24) 



we have again, as earlier, made the tacit assumption of neglecting the B"^ in comparison to the poloidal components. 
In the above equations, B — (Bp)^ + [B^Y is the magnitude of the magnetic field, gij and hi (with i,j = -0, p) 
are the metric coefficients and scale factors respectively and g is the determinant of the metric tensor. On a single 
field line, the metric reduces to: 



ds = hpdp 



(25) 



where s is the length along the magnetic field line. Equation (|25|) agrees with the equation of a magnetic field line 
(D'haeseleer et al. 1991). Integrating this last equation gives the length of the magnetic field line s,„ax(0'). 

In this new coordinate system, the azimuthal components of the Navier Stokes and magnetic induction equations 
reduce to: 



1 



h^BPhp 



d 



dt Aire ^Smax(V') ds 
dBf h^h„ d 



{h^B'P) 



dt 



{h^v'PBP) 



(26) 
(27) 



where we denote s 
below: 



[il)) ds ' 

s/smax(V'). These are the counterparts of the 2-D equations (in cylindrical coordinates) given 



dvf 
'df 

dBf 

~dr 



1 

'r 



1 

47re 
d 



d_ 

dp 



{v^BP 



(28) 
(29) 



2.1.2. Computational Domain 

The central singularity in the magnetic field now shifts to a coordinate singularity for ?/;. Numerically, we take care 
of this singularity by introducing a new variable q such that ij) — q/ (1 — q). We make another numerically convenient 
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Q. 




2 4 6 8 10 12 

s+(10xq) 



Fig. 1. The variation of p with respect to s + (10 x q) for a dipolar magnetic field. This and other figures (unless 
explicitly mentioned) employ 11 values of q; each value corresponding to a field line. For each line, the value of q can 
be read off from the x-axis of this figure (since 0<s<lasO<p< pmax, where p-max corresponds to the value of 
the cylindrical radial coordinate at the surface of the star). The lowest value of q corresponds to an 'open' field line 
almost parallel to the rotation axis while the largest one corresponds to a 'closed' field line enclosing the core of the 
star. The line corresponding to 5 = 0.5 is the marginally 'closed' field line intersecting the equator at the surface of 
the star. Predictably, the non-linearity of the variation increases with q. 



change of variable by integrating Eq. H25|l to obtain the length of the magnetic field line (smax) and then scaling 
the variable s with this quantity to obtain a new coordinate s = s/s^i^^. In Fig. ^ we display the result of such an 
integration and plot p with s + lOxq for the dipolar magnetic field. With these change of variables, the two independent 
spatial coordinates in the system of equations are q, and s such that < g < 1 and < s < 1. We have therefore, 
converted the curved field lines into straight field lines as shown in the left panel of Fig. [3 The projected field lines 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

s p 



Fig. 2. Computational domain: in the left panel are the field lines depicted by the 'crosses' and the 'pluses' for q 
v/s s; the vertical line represents the surface of the star (the region q > 0.5 represents the 'closed' field line that 
never intersect the stellar surface), the right panel displays the grid points projected onto the field lines in real 
space (cylindrical coordinates z v/s p) - the coordinate q labels each of these lines. The 'crosses' correspond to 'cell 
boundaries' at which B'^ is defined, while the 'pluses' correspond to 'cell centres' where ft is defined in the staggered 
space scheme that we use here. 
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with the grid points are shown in Fig. |21 right panel. Since the relationship between z and "0 is non-linear such that 
z blows up for low V' values, the central points for these values of ?/> look sparse in Fig. |21 right panel. 
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Fig. 3. Convergence tests: For the monopolar field case, for vacuum exterior. Top left panel is the convergence test 
for the minimum value of i?rot: the straight line corresponding to 0.513 is the analytically expected value; the limits 
on the y-axis represent 10% variation from this value, the '+' are the numerically obtained values with 51 grid points, 
the 'x' with 101 grid points, the '*' with 201, the filled squares with 401, the empty circles with 801 and the filled 
circles with 1601 (the last three appear merged in the main figure). The inset shows a blown up view of the same: solid 
line remains the same, the dashed line is with 101, the long dash-dot line with 201, the medium dash-dot line with 
401, the short dash-dot hue with 801 and the dot-dot-dashed line with 1601 grid points respectively. The top right 
panel is the maximum value of i?mag, the analytically expected value of 0.487 is represented by the straight line; the 
limits of the y-axis correspond to a 1% variation. The numerically obtained values for these are fairly well coverged 
even for 50 grid point resolution. The bottom left panel is the profile at a late time of = 29. The curves have the 
same meaning as for the i?rot convergence. The bottom right panel displays the B'^ profile sXt = 29. 



2.1.3. Initial and Boundary Conditions 

We assume that the differential rotation is zero at the center. This amounts to saying that: 



ds 



{t,q,s=Q) 



{t,q,s=Q) 



(30) 
(31) 



For the outer boundary conditions, let us first consider the 'closed' field lines. From symmetry about the equator, 
we know that all the derivatives in the vertical direction vanish in the equatorial plane, in addition = in the 
equatorial plane. From Eq. H29|) we see therefore that at s = 1 



B'^{t,q,s^ 1) = 



(32) 
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For the 'open' field lines, in a similar manner to Shapiro (2000), we can show using Roberts (1968) that for our 
straight field lines 



dt 



it,,,!) ds 



(l + 7^) „ 



(33) 

where TZ is the reflection coefficient for the azimuthal magnetic field amplitude: 



7^ = V^S^ (34) 

and VA = i?o/V47re is the surface Alfven speed with Bq = ^/[h^p]"^ being the representative surface magnetic field. It 
must be noted here that the actual internal Alfven speed, = B{p = 0, z = 0)/-\/47re >> i;a- 

For the initial conditions, as mentioned earlier, we assume that f2(0, q, s) is given by equation jnj. We also assume 
the toroidal magnetic fields to be absent to begin with i.e. B'^{Q, q, s) = 0. 

2.1.4. Non-Dimensional Formulation 

For the sake of numerical convenience, we scale the dimensions out of the equations through the following identifications: 

= [h^]h^, = [h^]h^, hp = [hp]hp, g = [g]g, v'^ = \h^^, 

n^QM, BP = B„BP, B^ = VI^enoK]B^, t ^ 2i/i[h^]/vA), IMMM = 

^ [hp] ) 

-^^It^ [hp] )[^^]^ 



2.1.5. Evolution Equations 

The non-dimensional form of the evolution Eqs. H2()|l and (|27|l become 
dn _ h^BPhp d 



di 2/i^V5Smax(g) 



{h^B'p) (35) 



dB'P h^hp d 



di 2^/gsm!ix{q) ds 



h^h^flBP) (36) 



where we have identified u*^ with h^pfl. 

2.1.6. Energy and Angular Momentum Conservation 

As energy conservation acts as a check to the numerical computations, we compute the relevant energies in the system. 
In the non-dimensional units, the integrals for rotational kinetic energy, toroidal magnetic energy and Poynting energy 
are given as (in the equations below, dr represents the volume element in cylindrical coordinates): 

£rot{t) - J dT{en{t,p,zfp^/2) = J d^E,,, (37) 

fmag(i) - / dT[{B^{t,p,z)f/%Tl]^ [ d^^E^^^ (38) 



£poy(t) = J dt£poy{t) = - J dt h^h^dil^dip ^ ^^^^^ ^ ^ J dt dipEpoy^ 

d^^Poy (39) 



S — Sri 
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where fpoy = —^h^B'-'^n and 

4ot = I f ds^g-^hl(f (40) 
J rip 

E^., = \ J ds^gk^iB-f (41) 



3 







^Poy = -2 / dt[B'PVLh^)^ (42) 



S — Sit 



are the energies associated with a given magnetic field line. Conservation of energy requires 

frot(O) = £rot(i) + fmag(t) + £poy (43) 

However, it may be noted that Eq. 1)43(1 is satisfied if on each magnetic field line: 

4ot (0) = 4ot (i) + ^mag {t) + E^oy {t) (44) 

holds. 

The angular momenta integrals are: 
J.ot[t) - j dT{en{t,p,z)p^) = J dxpl^t (45) 

J'mag(i) ^ J^dtAf = J d^A^mag (46) 



where A/" — £poy/ft is the torque exerted by the Maxwell stress at the surface and 



Jrot = ^ / dsv/^^/i^fl (47) 



P 

t 



Jmag = -J / di[B^h^) (48) 







2.2. Numerical Scheme 

Since the coefficients to the differentials in Eqs. (|35|) and (|35ll depend on some powers (> 1) of inverse of the coordinate 
distance from the centre, clearly, a fully explicit scheme may prove to be highly restritive due to the Courant condition. 
For the equivalent 1-D problem, it was shown (Shapiro 2000) that for a fully implicit scheme employing a 2nd order 
staggered space and 1st order time, the Courant condition only proves to be a restriction for accuracy (numerical 
stability being assured) - and the energy conservation (and angular momentum conservation in the case of energy 
dissipation due to an external plasma) is the marker for the accuracy. Accordingly, for our purpose here we use a 
staggered leapfrog scheme. The recipe for the finite differencing for such a scheme is mentioned in the Appendix. 

The details on the numerical checks performed on the code vis-d-vis the one-dimensional formalism of Shapiro 
(2000) is provided in section © below. 

3. Results 

In this paper, we calculate the timescales of the damping of differential rotation by magnetic fields. We have assumed 
a spherically symmetric star, having an applied (by hand) internal dipole magnetic field (that is divergence free 
everywhere and, in addition curl free everywhere except the centre). The initial differential rotation law is chosen 
(Shapiro 2000) to be given by Eq.© and initially, toroidal magnetic fields are assumed to be absent. For such a 
magnetic field configuration and initial rotation law, we solve the azimuthal components of the Navier-Stokes' equation 
and the Maxwell magnetic induction equation. For computational and physical convenience, we choose to solve these 
equations in a modified Magnetic Flux Co-ordinate system. This facilitates us to solve a series of 1-D equations, 
parametrised by a stream function. Using a staggered Leap-frog scheme, we difference the system of equations (as 
mentioned in Appendix) on a uniform q and s grid. In this section, we provide the results of our computations. For 
illustrative purposes, we have considered only 11 field lines (consistent with Figs. ^ and |2) and display results for 3 
field lines; in addition, we also display results for the surface and the equator (essentially one point for each field line) 
- for these, we have had to use 201 field lines for the purpose of clarity in the dynamics. 
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3.1. Checks on the code: simulations for a monopole held conhguration 

We first checked the code for a satisfactory reproduction of the 1-D results reported in Shapiro (2000). A trivial 
identification of {ip ^ z, ip ip, p p) gives us the monopolar magnetic field. In the above identification, the 
azimuthal Eqs. (|25|l and (|27|l reduce to those in Shapiro (2000) and as is mentioned therein, the energy conservation 
is an indicator of the accuracy of the numerical scheme. In Fig. |2| we display the convergence tests for the vacuum 
exterior case using resolutions corresponding to 51, 101, 201, 401, 801 and 1601 grid points on a single magnetic field 
line. The top left panel has the minimum values of E^ot (scaled to the value of the initial -Erot). The straight line located 
at 0.513 is the analytically obtained value for this quantity. The upper and lower limit on the y-axis correspond to a 
10% variation from this mean value. The top right panel displays the maximum value for the variations in E'mag (again 
scaled to the initial i?rot)- The analytical estimate of 0.487 is represented by the straight line located at the mid-point 
of the y-axis in the plot displaying a maximum variation of 1% in the limits. The bottom left panel displays the Q 
profile at time t = 29 a time where energy deviates substantially for the lowest resolution of 51 points. Understandably 
(from Eg. I40|l . the corresponding deviations in fl from the high resolution case is less significant than for i?rot- The 
bottom right panel displays the B"^ profile at t = 29. The ratio of the exterior density to the interior: Cox/e = 0.2. 

3.2. Simulations for the dipole field configuration 

We discuss the results for three field lines: one that is 'open' another that is 'closed' and yet another that is 'marginally 
closed' which we term equatorial, the latter two differing from the former through the boundary conditions as explained 
in section 1)2.1.3(1 and easily deduced from Fig. [3 In particular out of the 11 field lines that we consider here, we discuss 
the results for g = 0.3 (open), g = 0.6 (closed) and q = 0.8 (closed). 

3.2.1. Results for the open field line {q = 0.3) 

The dynamics of the variation of the rotation and the toroidal magnetic field are as shown in Fig.|31 We have used 851 
points along the field line direction. As can be seen, the torque between the stellar surface and the ambient plasma 
take away angular momentum from the star so that the rotation rate along the field line becomes uniform (Fig. [S] 
lower panel) . A plot of the energy and angular momentum variation over the dimensionless time t is depicted in Fig. 
Unlike for the monopolar case, the assumed ratio of density does not lead to retrograde rotation - this is because the 
curvature of the magnetic field line produces a lesser torque. 

3.2.2. Results for the closed field line {q = 0.6) 

The dynamics of the variation of the rotation and the toroidal magnetic field are as shown in Fig. |7| For this field 
line, we have used 1401 points in the line direction. Since this is a 'closed' field line and there is no dissipation in the 
system, the oscillations are sustained. At late times, the toroidal magnetic field profiles develop high gradients at the 
surface. This is due to d^l/ds not remaining zero at the surface at later times (as a result of the singularity of dz/dp 
at the equator). The energy conservation, however, still holds (the small variation from conserved value is due to the 
truncation error induced due to insufficient grid points). A key difference in the energy variation is the existence of 
a low frequency component (in comparison with the 1-D counterpart) - this difference is due to the fact that in the 
1~D case, the substitution of the induction equation into the Navier-Stokes equation yields a wave equation (Shapiro 
2000) in cylindrical coordinates, however, in our general form, such a substitution results in a wave equation with an 
advective component. The low frequency component in the energy variation with time, represents the oscillation of 
the point of inflection in our differential rotation law due to this advection. 

3.2.3. Results for the closed field line {q = 0.8) 

The dynamics (as earlier) is depicted in Fig. El The energy variation with time is shown in Fig. 1101 We have used 851 
data points. The deviation of the total energy for the value of 1 is due to truncation errors. 

3.3. Variation ofQ and at the surface 

We have seen from Figs. El [7| and that the dynamics of magnetic field and rotation rate variation with time along a 
given field line is smooth. However, if we use real (spherical/cylindrical) coordinates, the picture is different - since the 
oscillation along a field line is uncorrelated with that of an immediate neighbour, beyond a certain time, the oscillations 
get out of phase as can be seen in Fig. ^2 along the surface of the star. Structures form, however since Foynting fiux 
is carried away from the star, these settle down at late times to equilibrium values in off-equatorial regions. Since at 
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the equator, there are no toroidal magnetic fields generated, the oscillation is perpetual (unless diffusion or dissipation 
comes into play) - this fact, therefore, underlines the importance of taking into account diffusion to realistically model 
the system (Spruit 1999, Liu & Shapiro 2004). 

3.4. Variation ofQ in the equatorial plane 

Fig. 1121 illustrates the dynamical behaviour along the equator. The structures keep building up at late times here and 
it is important to consider diffusion or viscous dissipation so that the neighbouring lines 'talk' to each other. 

4. Discussion 

We have computed the dynamical evolution of rotation and toroidal magnetic fields for a spherical star composed 
of an incompressible homogeneous, differentially rotating fluid having an internal time independent dipolar magnetic 
field. Given these assumption one can expect the following physical processes to take place: At time t = 0, within 
the spherically symmetric differentially rotating star, a dipolar magnetic field is switched on. From Ferraro (1937) 
law of isorotation, we know that the lowest energy state of the system is when the rotation rate is constant along 
the magnetic field lines (provided there is dissipation in the system). For this to happen, the fiuid pulls the poloidal 
magnetic field lines, producing toroidal magnetic fields. In turn the toroidal magnetic field acts back on the plasma. 
In the absence of losses from the system, this leads to sustained oscillations. For an internal dipolar magnetic field, 
exchange between rotational and toroidal magnetic field energies will take place along the 'closed' field lines (fig) at 
characteristic frequencies. Along the open field lines, on the other hand, there will exist non-zero toroidal surface 
magnetic fields, that will be radiated away (or go into spinning up the ambient plasma) leading to losses from the 
system. This implies that after several Alfven time-scales, Ferraro law of isorotation will hold on the open field lines. 

We have shown here that for the 'open' field lines, rigid rotation is achieved (in accordance with Ferraro's law) 
through the loss of Poynting flux to an ambient plasma, while the 'closed' fleld line execute torsional oscillation through 
backreaction between toroidal magnetic and rotational kinetic energies. 

Although a first glance at Eq. [23 indicates that the right hand side becomes oo at the centre, it can be seen from a 
more careful observation and the acknowledgement of the fact that these equations are for a given magnetic field line, 
that the quantity within parenthesis will be and from our boundary conditions dfl/ds becomes zero at the centre. 
The upshot of this is that for different field lines, possesses different values at the centre - a manifestation of the 
singularity in the poloidal magnetic field and ip. 

An important result (and verification of the results of Liu & Shapiro 2004) is the relative difference in phases 
of oscillation of nearby points at 2-3 Alfven timescales when one samples a radial or polar zone (Figs. 1111 and [T^ : 
implying large phase deviation at late times (Spruit 1999). Such a behaviour is a result of our restrictive assumptions 
of neglecting diffusion and evolution of the poloidal field. Another point that the figures emphasise is the probable 
inefficacy of using pure spherical or cylindrical coordinates for analysing the dynamics of such a system in the absence 
of viscosity. 

The simplistic assumption of an internal dipolar magnetic field with 'closed' field lines is used as an illustrative 
example - the high restrictiveness of the sustenance of oscillation in such a system tests the efficacy of the code and the 
formalism in general. Hence, for the purposes of this paper, the origin and sustenance of such a field in matter is taken 
for granted. The formalism and the code amply demonstrate tenacity while the veracity of the numbers is ignored 
over obtaining the qualitative behaviour. The key aspect of the dynamics (when there is dissipation and according 
to Ferraro's law) is that it is all about matching two functions (through dissipation and readjustment of the angular 
momentum in the system): one that is constant along cylinders (the initial differential rotation rate) to another that 
follows the dipolar field. 

The evolution of the collapsing core of the progenitor star to pulsars is characterised by strong convection and 
turbulence. The convection is driven by neutrino transport and ceases when the star becomes transparent to neutrinos 
- an a-Lu dynamo operates in this regime (Thompson & Duncan 1993; Duncan & Thompson 1992) amplifying the 
primordial magnetic fields to saturation levels i3sat ~ 10"'^'* G. We dynamically place our system at a point after 
the cessation of convection and with a predominantly poloidal nature to the magnetic field - i.e. at about t^30 s 
after collapse. The next leg in the journey to computing equilibrium models of rotating magnetised stars will be the 
incorporation of less restrictive magnetic field geometries and the effects of dissipation. 
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Appendix A: Numerical Scheme 

In this section we describe the numerical scheme that we employ. The equations that we solve are of the type: 
dY d 

We assume X and Y to be functions of i, s and q (dependence on q is implicit as the eqns. ljA.l|) and ljA.2p are for 
constant values q). When we finite difference the above equations, we define X and Y on alternate grid points so as 
to yield a staggered mesh (and hence a second order treatment for space). Typically, we demand that X be defined 
on 'half grid points (cell boundaries) and Y on 'full' grid points (cell centres). Accordingly, the dependent variables 
will carry the tag: 



X(i,s,g) = Xr+i/^,,, Zr+3/2,, - e [1,^-1]) (A.3) 

Y{t,s,q) = Y^^, ... (^ e [\,N]) (A.4) 

where N is the total no. of grid points along a given magnetic field line (i.e. a given q). Since we're working on constant 
q (and hence constant j) lines, for the sake of easy representation here, we shall drop the subscript j from the variables 
(nevertheless remembering that the variables are indeed three dimensional). Thus, for a given q: 

x{t,s) = zr+3/2 ... e (A.5) 

Y{t,s) = Yl\ YIX^ ... {^ e [l,N]) (A.6) 
We discretise Eqs. (jA.l|) and (|A.2p in staggered leap-frog scheme as mentioned below: 

■"•^'--^'"^'^A.y.f "-"'^:"-^" ) (A,7) 



A.t \ As 



(A.8 
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Fig. 4. Convergence tests: For the monopolar field case, plasma exterior. In the top left panel is displayed the energy 

conservation, the right one displays angular momentum conservation while the bottom panel contains the CI and B'^ 
profiles at time t = 29. In all panels, the solid, the dashed, the short-dashed, the dot-dashed, the dot-short-dashed 
and the dotted lines indicate resolutions of 51, 101, 201, 401, 801 and 1601 points respectively. 
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Fig. 5. O and B"^ variation for the q = 0.3 case. The first panel shows the variation of Q as a function of the cylindrical 

coordinate p with time (the topmost curve being the initial curve and the rest the subsequent evolution). The second 
panel displays B'^ as a function of p with time (the initial and final profiles is a null profile and is hidden by the 
subsequent and earlier evolution respectively). The third panel shows the initial and final profile of Cl with p. The 
numbers labelling the curves in this and other figures stand for the corresponding time t; in this figure, since the curves 
get very close together at late times, we suppress labelling between t = 8 and t = 28. 
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Fig. 6. Energy and angular momentum conservation for the q — 0.3 case. The left panel displays the energy conserva- 
tion: the solid line is E^ot, the long-dashed line (coincident with the x-axis) represents £^mag, the short-dashed line is 
E'poy and the dot-dashed line is the sum of all these quantities (all energies are normalised to the initial value of Erot)- 
The right panel displays the angular momentum conservation: solid curve is Jrot (see text) while the short-dashed 
curve is the Jpoy and the dot-dashed line is the sum (all scaled to the initial Jrot)- 
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Fig. 7. n and B"^ variation for the q = 0.6 case. The curves have the same meaning as for Fig. |5l except that because 
q — 0.6 is a closed field line, the oscillations are sustained. 
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Fig. 8. Energy conservation for the q — 0.6 case. The curves have the same meaning as Fig.|Hl 
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Fig. 9. Vl and B'^ variation for the q — 0.8 case. The curves have the same meaning as for Fig.Oand dynamics almost 
identical to that for q — 0.6. 
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Fig. 11. il and B''^ variation at the surface. The curves have the same meaning as for Fig.[S| Structures form but do 
not build up due to the radiating away of the Poynting Flux. However, the dynamics at the equatorial surface (last 
point) do not cease due to this being a 'closed' field line. 
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Fig. 12. O variation at the equator; B*' is zero. Unlike the case for the surface, structures build up over time to a 
point where 'phase-mixing' become important. 



